function [X, L, X_err, L_err, iter, flag] = ...
    cheby(A, Xo, X_ref, L_ref, maxit, tol, print)
%
%  [X,L,X_err,L_err,iter,flag] = cheby(A, Xo, X_ref, L_ref, maxit, tol)
%
%  CHEBY applies the Chebyshev iteration to the eigenvalue problem
%
%    (A - L*I)*X = 0
%
%  The iterates take the form
%    X[l+1] = (1 / L[l]) * (A * X[l])
%    X[l+1] = X[l] + a[l+1] * (X[l+1] - X[l]) + b[l+1] * (X[l] - X[l-1])
%    L[l+1] = L[l] * ( A * X[l+1] ) / ( A * X[l] )
%  Reference:
%
%    Louis A. Hageman and David M. Yound
%    Applied Iterative Methods
%    Dover Publications, 2004.
%
%  Inputs:
%
%    real       A(N,N)     -- N-by-N matrix
%    real       xo(N)      -- Initial guess N-vector.
%    real       lambda_ref -- The actual eigenvalue for checking errors
%    integer    maxit      -- The maximum number of iterations.
%    real       tol        -- The error tolerance (defined...?)
%
%  Outputs:
%
%    real       x(N)       -- The computed eigenvector.
%    real       lambda     -- The computed eigenvalue.
%    real       rho        -- The estimated dominance ratio.
%    real       err_norm   -- The norm of the error.
%    integer    iter       -- The number of iterations performed.
%    integer    flag       -- flag (0=okay, 1=maxit reached, -1=bad bad news)
%

if (nargin==6)
    print = 1;
end

flag = 1;

% Reference values
% b   = -0.992856415680385/0.999999922699637;
% rho =  0.99403466159856/0.999999922699637;

x_L_minus_1 = Xo;
x_L         = Xo;
x_L_plus_1  = Xo;
lambda_L    = 1;

% cheb iterations using ref b and rho
p    = 0;
l    = 0;
iter = 0;
Npi  = 5;  % # power iterations for updating rho
tNpi = 20;  % total PI count
Nci  = 20; % # chebyshev iterations per cycle 
tNci = 0;  % total CI count
[x_L, lambda_L, rho] = powermethodRHO(A, x_L, tNpi); 

iter = tNpi;
while(iter <= maxit)
    
    if (p == 0)
        p = p + 1;
        iter = iter + Npi; tNpi = tNpi + Npi;
        [x_L, lambda_L, rho] = powermethodRHO(A, x_L, Npi);        
        b = -rho; % b is so close to rho in magnitude, this is probably fine
    end
    
    l = 0;
    while (l <= Nci && iter <= maxit)
        l = l + 1;
        iter = iter + 1;  tNci = tNci + 1;
        % update vectors
        x_hat_L_plus_1 = A * x_L / lambda_L;
%         rho = 0.993972526759680;
%         b   = -0.992776504518967;

        [alpha, beta] = alpha_beta(rho, b, l);
        
        x_L_plus_1     = (x_L + alpha*(x_hat_L_plus_1 - x_L) ...
                        +  beta*(x_L - x_L_minus_1)) ;
        
        % update eigenvalue        
        lambda_L_plus_1 = lambda_L*norm(A*x_L_plus_1,2)/norm(A*x_L,2);
        
        % save old 
        lambda_L    = lambda_L_plus_1;
        x_L_minus_1 = x_L;
        x_L         = x_L_plus_1;
        
        % I've gotta be able to effect this.
        [X, L] = powermethod(A, x_L, X_ref, L_ref, 2, tol, 0);
        
        X_err = norm(X - sign(X(1)/X_ref(1))*X_ref, 2);
        L_err = abs(L - L_ref);
        
        if ( X_err <= tol && L_err <= tol)
            flag = 0;
            break;
        end;
        if ( mod(iter, 50) == 0 )
            printout(iter,X_err,L,L_err,print,rho)
        end
    end
    if (flag == 0)
        break;
    end
    p = 0; % let another power loop occur
end

if (flag==0 && print == 1)
    disp(' *** Chebyshev: Final Results *** ');
end
if (flag==1 && print==1)
    disp(' *** Warning: Chebyshev did not converge ***');
end
disp(sprintf(' Power Iters = %5i, Chebyshev Iters = %5i, TOTAL Iters = %5i', ...
    tNpi, tNci, iter));
printout(iter,X_err,L,L_err,print,rho)
end

function printout(iter,X_err,L,L_err,print,rho)
if (print==1)
   disp(sprintf(...
       ' Iter = %5i, X_err = %6.5e, L = %6.4e, L_err = %6.3e, rho = %6.4e ',...
       iter, X_err,L,L_err,rho));
end
end

function [alpha, beta] = alpha_beta(rho, b, l)
% [alpha, beta] = alpha_beta(rho, b, l)
%
% Compute parameters for the Chebyshev recursion.  It is usually
% Here, b approximates lambda_N / lambda_1 and should be bounded from above
%  by that ratio.
% Here, rho approximates lambda_2 / lambda_1, the dominance ratio, and should be
%  bounded from below by that ratio.

% Reference values
%   b   = -0.992856415680385/0.999999922699637;
%   rho =  0.99403466159856/0.999999922699637;
if ( l == 1 )
    % special values for first iteration in cheby sequence
    alpha = 2 / (2 - rho - b);
    beta  = 0;
else
    y = (2.0 - rho - b) / (rho - b);
    alpha = 4.0 / (rho - b) * (Tn(l - 1, y) / Tn(l, y)); % Eqs. 11-2.12 in
    beta = Tn(l - 2, y) / Tn(l, y);                      % Hageman & Young
end

end

% Chebyshev polynomial of degree l evaluated at y
function z = Tn(l, y)
z = cosh(l * acosh(y));
end

% Estimate of rho
%   rho -- current estimate
%   Q   -- residual vector quotient
%   l   -- iteration index within current chebyshev cycle
function r = rho_est(rho, Q, l)
    r = 0.5*rho * ( cosh( acosh(Q*Tn(l,2.0/rho-1.0)) / l ) + 1.0 );
end

% Q-function from Hageman and Kellog.  Q is the residual vector
% quotient, and quantifies how much iterates are changing
function q = Q(R2,R1)
    q = norm(R2,1) / norm(R1,1);
end
